home *** CD-ROM | disk | FTP | other *** search
/ Amiga Tools 3 / Amiga Tools 3.iso / grafik / raytracing / rayshade-4.0.6.3 / fixes / fix024 / libray / libcommon / transform.c.diff < prev    next >
Encoding:
Text File  |  1994-08-09  |  5.0 KB  |  156 lines
  1. *** transform.c    Mon Oct 11 15:00:58 1993
  2. --- transform.c.frac    Mon Oct 11 14:47:18 1993
  3. ***************
  4. *** 305,310 ****
  5. --- 305,418 ----
  6.       MatrixCopy(&ttmp, inverse);
  7.   }
  8.   
  9. + /* find eigenvalues of positive definite matrices, return FALSE is singular. */
  10. + int MatrixEigenvalues(trans, lambda)
  11. + RSMatrix *trans;
  12. + Float lambda[3];
  13. + {
  14. +         Float a[3], maxval, a11, a12, a21, a22;
  15. +         Float *A[3];
  16. + 
  17. +         a[0] = fabs(trans->matrix[0][0]); maxval = a[0];
  18. +         a[1] = fabs(trans->matrix[1][0]); if (a[1] > maxval) maxval = a[1];
  19. +         a[2] = fabs(trans->matrix[2][0]); if (a[2] > maxval) maxval = a[2];
  20. +   
  21. +         if (maxval == a[0]) {
  22. +                 A[0] = trans->matrix[0];
  23. +                 A[1] = trans->matrix[1];
  24. +                 A[2] = trans->matrix[2];
  25. +         } else if (maxval == a[1]) {
  26. +                 A[0] = trans->matrix[1];
  27. +                 A[1] = trans->matrix[0];
  28. +                 A[2] = trans->matrix[2];
  29. +         } else {
  30. +                 A[0] = trans->matrix[2];
  31. +                 A[1] = trans->matrix[1];
  32. +                 A[2] = trans->matrix[0];
  33. +         }
  34. + 
  35. +         lambda[0] = A[0][0];
  36. +         if (lambda[0] == 0.) return FALSE;
  37. + 
  38. +         a[1] = A[1][0] / lambda[0];
  39. +         a[2] = A[2][0] / lambda[0];
  40. +         a11 = A[1][1] - a[1] * A[0][1];        
  41. +         a21 = A[2][1] - a[2] * A[0][1];        
  42. +         a12 = A[1][2] - a[1] * A[0][2];        
  43. +         a22 = A[2][2] - a[2] * A[0][2];        
  44. + 
  45. +         if (fabs(a11) > fabs(a21)) {
  46. +                 lambda[1] = a11;
  47. +                 if (lambda[1] == 0.) return FALSE;
  48. +                 lambda[2] = a22 - a21/a11 * a12;
  49. +         } else {
  50. +                 lambda[1] = a21;
  51. +                 if (lambda[1] == 0.) return FALSE;
  52. +                 lambda[2] = a12 - a11/a21 * a22;
  53. +         }
  54. + 
  55. +         if (lambda[2] == 0.) return FALSE;
  56. +         return TRUE;
  57. + }
  58. + 
  59. + /* Lipschitz number for euclidian norm: largest scaling factor */
  60. + Float Lipschitz(lambda)
  61. + Float lambda[3];
  62. + {
  63. +         Float maxlambda;
  64. + 
  65. +         maxlambda = lambda[0];
  66. +         if (lambda[1] > maxlambda) maxlambda = lambda[1];
  67. +         if (lambda[2] > maxlambda) maxlambda = lambda[2];
  68. +         return maxlambda;
  69. + }
  70. + 
  71. + /* returns TRUE if the three scaling factors are equal */
  72. + int Uniform(lambda)
  73. + Float lambda[3];
  74. + {
  75. +         return (equal(lambda[0], lambda[1]) && equal(lambda[1], lambda[2]));
  76. + }
  77. + 
  78. + /* scaling factors */
  79. + int MatrixScaling(trans, lambda)
  80. + RSMatrix *trans;
  81. + Float lambda[3];
  82. + {
  83. +         RSMatrix Q2;
  84. +         int i;
  85. + 
  86. + /* Q-R decomposition: we are only interested in Q**2 */
  87. +         Q2.matrix[0][0] = trans->matrix[0][0] * trans->matrix[0][0] +
  88. +                           trans->matrix[0][1] * trans->matrix[0][1] +
  89. +                           trans->matrix[0][2] * trans->matrix[0][2] ;
  90. +         Q2.matrix[0][1] = trans->matrix[0][0] * trans->matrix[1][0] +
  91. +                           trans->matrix[0][1] * trans->matrix[1][1] +
  92. +                           trans->matrix[0][2] * trans->matrix[1][2] ;
  93. +         Q2.matrix[0][2] = trans->matrix[0][0] * trans->matrix[2][0] +
  94. +                           trans->matrix[0][1] * trans->matrix[2][1] +
  95. +                           trans->matrix[0][2] * trans->matrix[2][2] ;
  96. +         Q2.matrix[1][1] = trans->matrix[1][0] * trans->matrix[1][0] +
  97. +                           trans->matrix[1][1] * trans->matrix[1][1] +
  98. +                           trans->matrix[1][2] * trans->matrix[1][2] ;
  99. +         Q2.matrix[1][2] = trans->matrix[1][0] * trans->matrix[2][0] +
  100. +                           trans->matrix[1][1] * trans->matrix[2][1] +
  101. +                           trans->matrix[1][2] * trans->matrix[2][2] ;
  102. +         Q2.matrix[2][2] = trans->matrix[2][0] * trans->matrix[2][0] +
  103. +                           trans->matrix[2][1] * trans->matrix[2][1] +
  104. +                           trans->matrix[2][2] * trans->matrix[2][2] ;
  105. +         Q2.matrix[1][0] = Q2.matrix[0][1];
  106. +         Q2.matrix[2][0] = Q2.matrix[0][2];
  107. +         Q2.matrix[2][1] = Q2.matrix[1][2];
  108. + 
  109. + /* scaling factors are square root of eigenvalues of Q**2 */
  110. +         if (!MatrixEigenvalues(&Q2, lambda)) return FALSE;
  111. +         lambda[0] = sqrt(lambda[0]);
  112. +         lambda[1] = sqrt(lambda[1]);
  113. +         lambda[2] = sqrt(lambda[2]);
  114. +         return TRUE;
  115. + }
  116. + 
  117.   /*
  118.    * Apply a transformation to a point (translation affects the point).
  119.    */
  120. ***************
  121. *** 399,407 ****
  122.   Ray *ray;
  123.   RSMatrix *trans;
  124.   {
  125.       PointTransform(&ray->pos, trans);
  126.       VecTransform(&ray->dir, trans);
  127. !     return VecNormalize(&ray->dir);
  128.   }
  129.   
  130.   void
  131. --- 507,519 ----
  132.   Ray *ray;
  133.   RSMatrix *trans;
  134.   {
  135. +         Float stretch;
  136.       PointTransform(&ray->pos, trans);
  137.       VecTransform(&ray->dir, trans);
  138. ! 
  139. !     stretch = VecNormalize(&ray->dir);
  140. !         ray->stretch *= stretch;
  141. !         return stretch;
  142.   }
  143.   
  144.   void
  145. ***************
  146. *** 433,437 ****
  147. --- 545,552 ----
  148.   {
  149.       TransCopy(list, result);
  150.       for (list = list->next; list; list = list->next)
  151. + /* it was: 
  152.           TransCompose(list, result, result);
  153. +  * PhB made it: */
  154. +         TransCompose(result, list, result);
  155.   }
  156.